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Abstract. We present Very Long Baseline Interferometry (VLBI) observations of the high mass X-ray binary 
LS I -|-61°303, carried out with the European VLBI Network (EVN). Over the 11 hour observing run, performed 
~10 days after a radio outburst, the radio source showed a constant flux density, which allowed sensitive imaging 
of the emission distribution. The structure in the map shows a clear extension to the southeast. Comparing our 
data with previous VLBI observations we interpret the extension as a collimated radio jet as found in several other 
X-ray binaries. Assuming that the structure is the result of an expansion that started at the onset of the outburst, 
we derive an apparent expansion velocity of 0.003 c, which, in the context of Doppler boosting, corresponds to 
an intrinsic velocity of at least 0.4 c for an ejection close to the line of sight. From the apparent velocity in all 
available epochs we are able to establish variations in the ejection angle which imply a precessing accretion disk. 
Finally we point out that LS I -|-61°303, like SS 433 and Cygnus X-1, shows evidence for an emission region almost 
orthogonal to the relativistic jet. 



> 

m 

a^ 
o 
l> 

o 

o 

:^: 

Oh 

O 

>- 1. Introduction 
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LS I H-61°303 is an X-ray binary system associated vi^ith 
the galactic plane variable rad io sou rce GT 0236-f 610 dis- 
covered by Gregory & Ta ylor (1978 ). Optical observations 
(Hutchings & Crampton 1981 ) show that the system is 
composed of a neutron star and an early type, rapidly ro- 
tating BOV star with a stable equatorial disk and mass 
loss. Spectral line observations of the radio source give a 
distance of 2.0±0.2 kpc (Frail & Hjellming |l99ll ). 

One of the most unusual aspects of its radio emission 
is the fact that it exhibits two periodicities: a 26. 5 day 
perio dic nonthermal outburst (Taylor 

4 years 
al 



Simultaneous X-ray and radio observations show that 
the X-ray outbursts occur at the periastron passage while, 
on the contrary, the strongest radio outbursts are always 



delaye d wit h respect to that (Taylor et al. 1996; Gregory 
et al. 1999| ). Both, the presence of two periodicities (at 
26.5 days and ~ 4 years) and the delay between radio 
and X-ray outbursts, are well explained in the framework 
of an accretion scenario of a precessing neutron s tar in 
a highly (e > 0.4) eccentric orbit (Gregory et al. 1989 ; 



Gregory et al. 1999). The accretion rate in an eccentric 



1982 



198J) and a 1584 day (- 



Gregory 
mo dulation of the out- 
burst peak flux (Gregory et al. |1999 ). The 26.5 day pe- 
riodicity corresponds to the orbit al pe riod of the binary 
system (Hutchings & Crampton 1981). This periodicity 
has also been detected in UBVR I photometric observa- 
tions (Mendelson fc Mazeh |1989|) , in the infrared domain 



(Paredes et al. |1994| ), in soft X-rays (Paredes et al. |1997D 
and in the Ha emission line (Zamanov et al. 1999). The 
4 year modulation has been observed as well in the Ha 
emission line (Zamanov et al. 1999). 
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orbit within the equatorial wind of the Be star has two 
peaks: the highest peak corresponds to the periastron pas- 
sage and the second, lower amplitude, peak occurs when 
the relative velocity of the neutron star and the Be star 
wind is at a minimum. For supercritical accretion, mat- 
ter is ejected outwards in two jets perpendicularly to the 
accretion disk plane. Near periastron, inverse Compton 
losses are severe (due to the proximity to the Be star): 
X-ray outbursts are expected but not radio ones. For the 
second accretion peak, the neutron star is much farther 
from the Be star and both inverse Compton losses and 
wind opacity are lower, the electrons can propagate out of 
the orbital plane and we observe the radio outburst. The 
precession of the disk gives rise to the ~ 4 year modulation 



(Gregory et al. 1989; Taylor et al. 1992: Massi et al. 1993 
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Fig. 1. Radio light curve of LS I +61 °303 o btained with 
the GBI at 8.4 GHz (Strickman et al. |1998D . The shaded 
area indicates the time interval during which our VLBI 
data were taken 



Marti & Paredes 1995; Gregory et al. 1999). The pres 



ence of an accretion disk is also invoked by Mendelson & 



Mazeh (1989) to explain details of the optical light curve. 
Liu et al. (200C) explain the variation of the Ha emission 
with orbital phase as varying irradiation of the Be star's 
circumstellar disk by the X-ray emission from the neutron 
star's accretion disk. 

However, the luminosity of LS I +61°303 in the X-ray 
range ( 1-40 keV), is only Lx — 10'^^ erg s^^ (Maraschi & 



Treves 1981 ), which is three orders of magnitude lower 
than the Eddington limit. Moreover, LS I +61°303 is 
the most promising candidate for the optical counterpart 
of th e 7-ray source BEG J0241+6103 (Hartman et al. 
1999| ), with a luminosity L^ ~ 10"^^erg s~^. The fact that 



LS I -|-61°303 has the bulk of its energy shifted from X- 
ray to 7-ray wavelengths is not understood in the context 
of the supercritical accretion model, showing the neces- 
sity for alternative models. These models suggest that the 
26.5 d outburst events are produced by energetic electrons 
accelerated in the shock boundary between the relativis- 
tic wind of a young non-accreting pulsar and the wind of 
the Be star. In this case the 4 year modulation would be 
explained due to cyclic variations of the Be star envelope 
(Maraschi & Treves |1981| Ta vani |1994 Tava ni et al. |1998 ; 
Goldoni & Mereghetti |1995| ; Zamanov [1995| Taylor et al. 
199^ ). 

The recent discovery of the microquasar LS 5039 
(Paredes et al. 2000) brings new credibility to the ac- 
cretion model. LS 5039 is also subluminous in the X-ray 
range (even more than LS I -|-61°303) and also shows 
the same puzzling behavior, having Lj > Lx- Therefore, 
LS 5039 and LS I +61°303 could be the first two exam- 
ples of a new class of X-ray binaries with powerful 7-ray 
emission. In the case of LS 5039, an ejection process fed 
by an accretion disk is clearly proved by a map obtained 
with the Very Long Baseline Array (VLB A), which shows 
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Fig. 2. Effelsberg flux density measurements of 
LS I +61°303 at 5 GHz obtained during our VLBI 
observations. LCP and RCP flux density measurements 
are represented by filled and open circles respectively. 
The solid line corresponds to the total flux density 



bipolar jets emerging from a central core. Its high L-y is 
tentatively explained by inverse Compton scattering. In 
LS I +61°303, although several VLBI observations show 
a complex source extending over a few milliarcseconds 
(Massi et al. 



1993; Peracaula et al. 199 J; Paredes et al. 



|1998 ; Taylor et al. 2000), such a clear jet structure has 



never been observed. 

In this paper we report on EVN observations probing 
structures on scales of tens of milliarcseconds, larger than 
probed by previous observations. We clearly detect on this 
scale an elongation of the emission (Sect. El), which we in- 
terpret in Sect. 0, as a one-sided jet. We determine its 
apparent expansion velocity and derive the intrinsic ve- 
locity considering Doppler boosting. Our conclusions are 
given in Sect. 0. 



2. Observations and results 

Our observations of LS I +61° 303 were carried out during 

II hours on 1994 June 7, at 5 GHz, using four anten- 
nas of the EVN, namely the Effelsberg 100m, Medicina 
32m, Noto 32m and Onsala 25m telescopes. The data were 
recorded, in left hand circular polarization, using the Mark 

III mode A recording system, corresponding to a total 
bandwidth of 56 MHz. The correlation of the data was 
done on the Mark HI correlator operated by the Max- 
Planck-Institut fiir Radioastronomie in Bonn. 

Our EVN observations were made within a time inter- 
val during which a multiwavelength monitoring campaign 
was conducted that included simultaneous radio, optical, 
infrared, and hard X-ray/7-ray observations and covered 
almos t thre e orbital cycles of LS I +61°303 (Strickman 
et al. 1998|). As can be seen from the radio light curve 
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taken with the Green Bank Interferometem (GBI) as part 
of that campaign (Fig. nl), our observations took place 
about 10 days after the radio outburst maximum (i.e. 
almost 13 days after the onset of the outburst). Total 
flux density observations of LS I +61°303, carried out by 
the 100m telescope of Effelsberg at the beginning of each 
VLBI scan, are shown in Fig. 0. The variations of left and 
right hand circular polarization flux are within their error 
bars. This lack of significant circular polarization agrees 
with previous results of Peracaula et al. (1997). The total 
flux density is nearly constant and with an average value 
of 34 ± 5 mJy. 

In our EVN observations, LS I +61°303 was clearly 
detected on all baselines. The data were reduced using 
standard procedures within the AIPS and Difmap soft- 
ware packages. The EVN uniform weighted map is shown 
at the top of Fig. ^. It exhibits an extended structure 
elongated to the southeast. The elongation is not artifi- 
cially created by the beam, in fact the elongation is at 
P.A.~ 120° while the beam, of size 5.9 mas x 3.8 mas, 
has P.A.=74.2°. The flux density recovered in the clean- 
ing process amounts to ~ 35 mJy. 

In Table |l| we show the parameters obtained after fit- 
ting several models to the (u, v) data. The best fit using a 
single component is obtained with an elliptical Gaussian 
(Model 1 in Table |l|) . A circular Gaussian or a point source 
give higher values of Xr- The elliptical Gaussian has a 
P.A. consistent with the observable elongation in the map. 
However, there is still residual emission not accounted by 
Model 1, in the direction of the elongation. 



In previous VLBI observations (Paredes et al. 1998) 
the best fits to the (u, v) data were consistent with an un- 
resolved core plus a Gaussian halo. Hence, we have added 
a point source to Model 1 , and we have obtained a slightly 
improved fit to the data (Model 2 in Table mj. It is im- 
portant to note that, in this model, the elliptical Gaussian 
conserves its P.A. and moves towards the direction of the 
elongation, while the point source moves in the opposite 
direction. We have tried several different initial conditions 
and always the fit converged to the same soluti on. Fi nally, 
for a distance of 2.0 kpc (Frail & Hjellming 1991 ), the 
brightness temperature of the components of the models 
listed in Table |]are in the range of Tb = 10^-10^ K, which 
are characteristic of non-thermal emission. 

3. Discussion 



3.1. Morphology 

In view of the models that can fit the data, the source in 
Fig. y can be explained as a central core (point source of 
Model 2 in Table |l]) and a one-sided jet (elliptical Gaussian 
in the same model). If we assume that the Gaussian com- 
ponent is the result of an expanding source that originated 
at the onset of the major outburst, that is about 13 days 

^ The Green Bank Interferometer is a facility of the USA 
National Science Foundation operated by the NRAO in sup- 
port of NASA High Energy Astrophysics programs. 
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Fig. 3. Top: EVN uniform weighted map of LS I +61°303 
at 5 GHz frequency obtained on 1994 June 7 (this paper). 
The contours are at -3, 3, 4, 6, 8, 11, 16, 22, 30, 45, 60 and 
75 times the r.m.s. noise of 0.28 mJy beam~^. The filled 
ellipse in the bottom-left corner represents the FWHM 
of the synthesized beam, which is 5.9 mas x 3.8 mas at a 
P.A. o f 74.2 °. Bottom: VLBI map (at 5 GHz) by Massi 
et al. (1993), obtained on 1990 June 6, with a resolution 



of 0.6 mas x 0.5 mas 



before our observations, and that the expansion velocity 
of the jet is constant, the derived velocity, taking into ac- 
count the size of the Gaussian and the offset to the point 
source, is ~ 0.3 mas d~^. This corresponds, at a 2.0 kpc 
distance, to a projected expansion velocity of ~ 0.003 c. 

The most interesting aspect of the EVN map is that, 
for the first time, we detected asymmetric emission in the 
southeast direction. The morphology is quite similar to 
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Table 1. Parameters of different models fitted to the data 



Model 


xl 


Components 


Flux Density 
(mJy) 


r 

(mas) 


e 

n 


Maj. axis 
(mas) 


Min. axis 
(mas) 


P.A. 

n 


1 


0.798 


Elliptical Gaussian 


33.0 


0.04 


3 


3.7 


2.5 


123 


2 


0.786 


Elliptical Gaussian 
Point source 


21.6 
13.1 


0.5 
0.5 


135 
-43 


5.5 


3.6 


123 



Table 2. Summary of published VLBI observations of LS I +61°303 at 5 GHz. Phase is computed according to 
P = 26.4917 d and Tq = JD 2 443 366.775 



Epoch 


Array 


Phase 


State 


Flux Density 
(mJy) 


Size 
(mas) 


P.A. 

(°) 


Expansion 
Velocity/c 


1987 Sep 25" 
1987 Oct 1" 


EVN 


0.59 
0.81 


quiescent 
burst, decaying 


40-54 
260-200 


3.2 ±0.9 
1.6 ±1.2 




< 0.002 


1990 Jun 6'' 


EVN + VLA 


0.74 


burst, decaying 


244-205 


~2 


~ 135, ~ 30 


~ 0.002^ 


1992 Jun 8" 


Global VLBI 


0.42 


quiescent, minioutburst 


35 


0.5-2 


~ 160 


0.06 


1993 Sep 9^* 
1993 Sep 13'' 


Global VLBI 


0.69 
0.84 


burst, variable 
quiescent 


76-131 
60 


2-3 

~7 


~ 120 


0.007 


1994 Jun 7" 


EVN 


0.92 


quiescent 


34 


rsj 6 


~ 120 


0.003 


1999 Sep 16-17^' 


HALCA + 
Global VLBI 


0.69 


burst, variable 


140 


rvj 4 




~ 0.002^ 



° Taylor et al. |l992| ; '' Massi et al. [l993| ; "= Peracaula et al. [l998| ; "^ Paredes e t al. J l998| ; " This paper; ^ Taylor et al. |200C 



^ Recalculated according to new outburst peak ephemerides (Gregory et al 

^ There is not enough resolution at P.A. ~ 140 to know if there was expansion in that direction 



t al.F 
1999^ ^ 



1999), for the extended structure at P. A. =135° 



that observed by Stirhng et al. ( pOOO , Fig. 4) in Cygnus X- 
1 at 15 GHz: a central source with a small but clear elon- 
gation. Is this morphological analogy enough to interpret 
the structure in the EVN map as a one-sided jet? In the 
case of Cygnus X-1 the interpretation of the small jet-like 
elongation has proved to be correct: a later observation of 
this source at 8 GHz shows the elongation developed in 
an extended jet. 

LS I +61°303 has been observed several times with 
mas resolution. In Table we report previous VLBI ob- 
servations focusing on two special items: the position angle 
(P.A.) of any extended feature present and its expansion 
velocity. As can be seen, an extended feature with P.A. 
from 120° to 160° seems to be present in other maps. 
While the identification is perhaps ambiguous in Paredes 
et al. 1998 because the extension is coincident with the 



direction of the major axis of the beam, the observations 
by Peracaula et al. (1998), report on one extended fea- 
ture with comparable P.A. with that of the feature in our 
map. The reason why in the past this extended structure 
was never clearly associated to a jet is because an internal 
structure at another P.A. (~ 30°) makes the morphol- 
ogy confusing. To better understand this point we show 



at the bottom of Fig. the map by Massi et al. (1993). 
The VLBI map resolves two components (at P.A.~30°) in- 
side the extended structure (at P.A.~ 135°). The extended 
structure is asymmetric as well as that in our map, but 
towards the northwest direction. This is probably a pro- 
jection effect corresponding to a different ejec tion a ngle 
due to precession (Sect. ^^). In Paredes et al. ( 199^ ) the 
double internal source is not resolved, but the authors pro- 
pose the presence of an unresolved component stable and 
not participating in the flaring process, inside an extended 
structure at P.A.~ 120°. These two components are again 
resolved in the recent VLBI observation using HALCA 
(Taylor et al. 2000). While in the past this double inter- 
nal source was supposed to be the jet, in light of the EVN 
map we propose that: 



The double source at P.A. of ~ 30° is not the jet. 
Instead, it is almost orthogonal to the jet. It has a size 
of 2 AU and therefore is comparable with the orbital 
size (1.4 AU). 

The jet is identified with the structure with P.A. in the 
range 120°-160° seen in at least three VLBI observa- 
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tions, and best defined in our EVN map where it has 
a size of ~ 8 AU. 



An emitting region perpendicular to the jet has been 
observed indeed also in SS 433 by Paragi et al. ( 1999 ) 
and perhaps in Cygnus X-1 by Stirhng et al. ( 2000 ) 
and preliminary interpreted as shocked gas in the or- 
bital/accretion plane. If this interpretation holds true also 
for LS I +61°303 the strong wind opacity close to the pe- 
riastron, discussed in the introduction, certainly severely 
affects the morphology of such emission. For the jet, how- 
ever, wind opacity effects are negligible. The jet is gener- 
ated by outbursts occurring at quite a displaced orbital 
phase with respect to the periastron. Moreover, the size 
of the jet is much greater than the orbital size. 



3.2. One-sided Jet and Doppler boosting 

It is well known that while some extragalactic radio 
sources have two jets, for some others only one jet is ob- 
served. The unification model for AGN (see a review in 
Urry & Padovani 1995) assumes that all of them repre- 
sent the same class of objects (all having two jets) and 
their different appearance depends on different observing 
angles. Let us assume a symmetric ejection of two jets at 
velocity /3 (i.e. expressed as fraction of c). The two jets, 
approaching and receding, move at an apparent velocity 
/3a,r related to the intrinsic /3 by (Rees 1966; Mirabel & 
Rodriguez [1994D 

_ /3sing 
^-■^-l^/3cos0' ^' 

where 9 is the angle between the direction of motion of 
the ejecta and the line of sight. 



Following the method of Mirabel & Rodriguez (1994) 
one can determine the quantity /3 cos 6 by means of the 
ratio of flux densities from the approaching and receding 
jets. 



^ 

^r 



1 + /3 cos 6* 
I- f3 cos 



(2) 



where a is the spectral index of the emission (S oa i>°') and 
fe is 2 for a continuous jet and 3 for discrete condensations. 
In our case we deal with one jet only. However, we can 
determine the lower limit 



13 COS 9 > 



1 



{Sr^^ySa 



l/{k-a) 



1 



(3) 



using the noise level (a) of the map and the peak value of 
the approaching component {Sl°^^). The spectral index a 
few days after the outburst is a = —0.5 (Strickman et al. 



199a). To be consistent with the lowest limit we select 



fc = 3. 

For our EVN map, with a = 0.28 mJy beam^^ and 
^poak ^ 21.6 mJy beam~\ we have (3cos9> 0.43. 
Using Eq. ^ we obtain 



which for /3a = 0.003 (Sect. |3T|) and (3 cos 9 > 0.43 gives 
an intrinsic velocity of /3 > 0.43 and an ejection angle 
6* ~ 0°. If we consider a typical size of ~ 10~^ pc for the 
jets in other known microquasars, the ~ 8 AU size of the 
LS I -f61°303 jet would result in an angle of 6^ ~ 0?2, 
compatible with our estimation. 

As reviewed by Mirabel & Rodriguez ( |1999D , the ex- 
pansion velocities for microquasars range from ^ 0.1 
c to - 0.9c (SS 433, Cygnus X-3, GRS 1915-fl05, 
GRO J1655— 40). The above determined lower limit of 
0.4c for LS I -|-61°303 is therefore well within that range. 
Finally, we note that the combination of the values esti- 
mated for 9 and (3 gives a Lorentz factor 7 — 1.1, and a 

r, (5a — 1-59 for the ap- 



Doppler factor 5^^, = ^(i^^co.e) - 

proaching jet and 5^ = 0.63 for the receding one. 



3.3. Precession of the accretion disl< 

The measured expansion velocities shown in Table |^ span 
a range of 0.002-0.007 c and reach a value of 0.06 c at 
epoch 1992 June 8. On the basis of our discussion in 
Sect. ^T] and 3^, we interpret the observed expansion ve- 
locities as apparent transverse velocities, /3a, defined by 



(3^ \l{(3cos9f +(3^^ {I- p cos 9f, 



(4) 



Eq. y. A possible explanation for the large range of ob- 
served velocities is a variable intrinsic velocity /3 and a 
constant 9. However, this is not supported by the obser- 
vations available up to now; moreover, for example, SS 433 
has shown a quite constant velocity of /3 = 0.26 for years. 
The alternative explanation is precession of the jet, with 
the angle 9 between the direction of the jet and the line of 
sight being a function of time. Evidence for precession has 
been found at least for SS 433 and Cygnus X-1 (Brocksopp 
et al. 1999| ). Moreover, for LS I -|-61°303, precession of the 
jet has already been suggested to explain the 4 yr modu- 
lation of the peak of the radio outbursts (Gregory et al. 
19891) . 

If the latter assumption is correct, we would expect 
an anticorrelation between the flux density of the radio 
outburst peak, and 9 or /3a. In other words, when the 
jet is pointing directly towards us {9 small), /3a is also 
small, and the flux density of the outburst peak is the 
highest possible due to the Doppler boosting effect. On 
the contrary, when the jet is not pointing directly towards 
us, 9 and /3a increase, and the flux density decreases. In 
Fig U we show the flux density of the radio outburst peak 
versus log /3a, for the VLBI observations listed in Table |[ 
taking /3a as the observed expansion velocity. We can see 
from Fig Q that the available data are consistent with the 
anticorrelation of the flux density of the radio outburst 
peak and /3a predicted by our model, and thus with the 
precession of the jet of LS I +61°303. 

Assuming that the the intrinsic velocity of the jet re- 
mains constant (/3 ~ 0.4) for all the observations listed 
in Table g, we can derive that the maximum value of the 
ejection angle 0, corresponding to the maximum value of 
/3a = 0.06, is 9 ~ 4? 5. This maximum angle is consistent 
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4. Conclusions 




Fig. 4. Flux density of the radio outburst peak versus 
log /3a (see text). Each data point in this diagram is the 
result of the corresponding VLBI observation identified 
with letters a to f (see Table |]). For observations a, b, 
e and / the flux density of the radio outburst peak was 
known from the same VLBI observations or from the GBI 
monitoring program. Point d is lower limit of the outburst 
peak, while c is taken from the average outburst peak flux 
modulation of Gregory et al. (199£) 



Our EVN observations at 5 GHz revealed that 
LS I +61°303 has a jet-like elongation to the southeast 
on scales of tens of milliarcseconds. This is interpreted 
as an approaching, Doppler boosted, jet of a symmetric 
pair. The angle between the jet and the line of sight at 
the epoch of our observation was close to zero. However, 
data at other epochs suggest the existence of precession 
of the jet, and the accretion disk, of at least 5°. This pre- 
cession would change the angle of the ejecta and conse- 
quently would vary both the apparent expansion velocity 
and the intensity of the peak of the radio outbursts. New 
VLBI observations at several epochs are necessary to cor- 
relate morphology, position angle, and expansion velocity, 
with the 4 yr modulation of the peak of the radio out- 



bursts (Gregory et al. 1999) and of the Ha emission line 
(Zamanov et al. 1999). 

We derived a lower limit of 0.4 c for the intrinsic 
velocity of the radio jet. This value is well within the 
range 0.1 c to 0.9 c found for microquasars like SS 433, 
Cygnus X-3, GRS 1 915-H05 and GRO J1655-40 (Mirabel 
& Rodriguez [l999| ). 

Finally, in some microquasars emission almost orthog- 
onal to the jet, i.e. along the orbital plane, exists. This 
emiss ion has clearly been observed in SS 433 (Paragi et al. 



1999) and it is barely visible in Cygnus X-1 (Stirling et al. 



2000). As we showed in this paper, LS I -|-61°303 could be 



third case of this kind of sources. 



with a moderate precession of the jet, and could be the 
result of the precession of the accretion disk. 

The range of values obtained for 9 imply that the 
orbital plane of the system would be close to the plane 
of the sky, i.e., the inclination of the orbit, i, should be 
small. The other available information about i comes from 
spectroscopic optical and UV observations carried out by 



Hutchings & Crampton (1981). Although their observa- 
tions reveal shell absorption, this fact alone does not give 
information on the angle of inclination, because the shell 
can cover the whole star and n ot be just an equatorial 
bulge (Kogure 1969| ; Geuverink 1970| ). We note that the 
observations by Waters et al. ( |l988|) establish that, to- 
gether with a dense and slow disk-like wind around the 
equator, there exists a high velocity, low density wind at 
higher latitudes up to the polar regions. On the other 
hand, Hutchings & Crampton (1981) obtain wsin i = 
360 ± 25 km s~^, where v is the equatorial rotational ve- 
locity of the Be star. As the maximum rotational velocity 
of a Be star (Hutchings et al. 1979 ) is w = 630 km s"^, 
this would result in a value of i ~ 35°. However, Hutchings 
& Crampton ( |l98l| ) comment that the velocity data are 
" extensive and unwieldy" and indeed the value of v sin i 
given above seems large if compared with the statistical 
study of Be stars by Slettebak (1982). New observations 
would be very useful to clarify this issue. 
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